%==========================================================================
% Matlab code to replicate the results discussed in online Appendix F.5 of 
% "Marital matching, economies of scale and intrahousehold allocations" 
% by Laurens Cherchye, Bram De Rock, Khushboo Surana and Frederic Vermeulen 
%==========================================================================

%% loading the data
clear all
load('psid_cleaned_withsingles.mat')
leisurem(isnan(leisurem))=0;
leisuref(isnan(leisuref))=0;
hworkm(isnan(hworkm))=0;
hworkf(isnan(hworkf))=0;
wagem(isnan(wagem))=0;
wagef(isnan(wagef))=0;

% adding paths
% replace the directory with the path where YALMIP is downloaded
addpath(genpath('C:\Users\u0100562\Downloads\YALMIP-master'))

%% Bounds on economies of scale and RICEB measures

clc
ub = 0.0;
bounds = [];
temp_id = id;

while isempty(temp_id)==0  
    fprintf('Working on marriage market id %d \n', temp_id(1));
    %preparing the id index for every id
    id_ind = find(id==temp_id(1));
    display('Computing divorce costs');
    [temp_DivorceCostsM,temp_DivorceCostsF,temp_DivorceCostsMF] = get_divorce_costs(singlefemale(id_ind),singlemale(id_ind),couple(id_ind),leisurem(id_ind),leisuref(id_ind), hworkm(id_ind), hworkf(id_ind), q(id_ind),wagem(id_ind),wagef(id_ind),nonlabor(id_ind),ub);
    display('Computing bounds');
    [temp_min_scale,temp_max_scale,temp_min_scalem,temp_max_scalem,temp_min_scalef,temp_max_scalef] = get_scale(singlefemale(id_ind),singlemale(id_ind),couple(id_ind),leisurem(id_ind),leisuref(id_ind), hworkm(id_ind), hworkf(id_ind), q(id_ind),wagem(id_ind),wagef(id_ind),nonlabor(id_ind),temp_DivorceCostsM,temp_DivorceCostsF,temp_DivorceCostsMF,ub);
    
    % merging bounds together 
    bounds_temp = zeros(length(find(couple(id_ind))),11);
    bounds_temp(:,1) = temp_min_scale';
    bounds_temp(:,2) = temp_max_scale';
    bounds_temp(:,3) = temp_max_scale' - temp_min_scale';
    bounds_temp(:,4) = temp_min_scalem';
    bounds_temp(:,5) = temp_max_scalem';
    bounds_temp(:,6) = temp_max_scalem' - temp_min_scalem';
    bounds_temp(:,7) = temp_min_scalef';
    bounds_temp(:,8) = temp_max_scalef';
    bounds_temp(:,9) = temp_max_scalef' - temp_min_scalef';
    bounds_temp(:,10) = familyid13(id_ind(find(couple(id_ind))));
    bounds_temp(:,11) = ub*ones(length(find(couple(id_ind))),1);
    bounds = [bounds; bounds_temp];
    
    temp_id(find(temp_id==temp_id(1)))=[];
    
end

fprintf('Computations for ub equal to %d are finished \n', ub);
save psid_results_f5_withsingles_different;




